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1. Introduction 


Ionospheric scintillations are fluctuations in amplitude and phase of the radio wave signal 
passing through the ionosphere on its path to the receiver, caused by small scale irregularities 
in the electron density structure [1]. Occurrence of scintillation depends on various factors, 
such as: solar activity, geomagnetic conditions, geo-location of the receiver, time of day, angle 
of signal arrival, ionospheric structure geometry and many others [1-6]. Causing signals 
disturbances and distortion, scintillation can significantly affect the GNSS accuracy and cause 
severe problems to commercial navigation systems. Scintillating signals can be classified by 
the intensity of fluctuations into categories of weak and strong scintillations. High and low 
magnetic latitudes represent the most affected regions by ionospheric irregular structures in 
F and E ionospheric layers. While phase scintillations are more pronounced in the sub-polar 
and polar regions, amplitude scintillations are significantly stronger and more pronounced in 
near equatorial regions [7-10]. 


In last five decades, statistical studies of the ionospheric scintillating radio signals brought 
different approaches and solutions. Some of these solutions are widely adopted and used in 
scintillation modelling and forecasting ionospheric dynamics, as joint Gaussian distribution 
of complex radio wave signal. This solution gives applicable results for weak scintillating 
signals [3], but in case of strong scintillations there is not an easy way to derive satisfactory 
results, leading to a need for further investigations [8]. While phase scintillations mainly follow 
the Gaussian distribution and do not represent a problem in ionospheric modelling, more 
complicated case is with amplitude scintillations. Revolution in ionospheric modelling and 
developing of the ionospheric scintillation theory has been made during '70s and ‘80s, covered 
mostly in papers by Rino and Fremouw [11-16]. During this period various researches had 
been performed providing different information on a probability distribution function (PDF) 
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used as amplitude scintillation descriptor [11-12, 17]. Last two decades brought few intriguing 
studies [18-21] performed on statistics of scintillating signal leading to aneed for more detailed 
and precise description of radio scintillation signal’s amplitude and phase PDFs for strong 
fluctuations. 


The focus of this paper is on the statistical analysis of ionospheric scintillation in high and low 
latitude during strong and moderate geomagnetic activity. Gaussian and Nakagami-m PDFs 
of scintillating signals have been examined using real measured data and compared with 
theoretically derived PDFs. Further testing had been done on different data intervals, neces- 
sary for correct higher order statistical analysis, avoiding errors influence. The analysis results 
are presented with higher order moments, dependent on various parameters (scintillation 
indices, geo-location and solar/magnetic activity). Implementation of higher order moments, 
skewness and kurtosis, could give additional information about the ionospheric irregularities 
influence on the propagating signal and relation to the time delay of the signal. 


2. Analysed data and methods 


The data used in analysis have been measured by NovAtel GPS lonospheric Scintillation/TEC 
Monitor (GISTM), model GSV4004B, one of three monitoring receivers installed at the Polish 
Polar Station in the region of the Hornsund Bay in southern Spitsbergen (approximately 77° 
N, 15.55° E). The second data source was derived using a Septentrio's PolaRxS Ionospheric 
Scintillation Monitoring (ISM) receiver located at Presidente Prudente (approximately-21.99° 
N, 308.59° E), Brazil, set during Concept for lonospheric Scintillation Mitigation for Professio- 
nal GNSS in Latin America (CIGALA) project. The geographical locations of both stations are 
indicated in the Fig. 1. The locations of the receivers made possible study of the ionospheric 
scintillations in the most affected regions characterized by very intense and complex behaviour 
of plasma - polar and auroral latitudes, due to interaction of the solar wind and interplanetary 
magnetic field with the Earth's magnetic field, and low latitudes, due to equatorial anomaly 
and electrojets. Parameters of interest in performed data analysis are signal amplitude and 
phase during quiet and active geomagnetic periods, recorded during measurement cam- 


paigns. 


Both receivers collect raw amplitude and phase output data with a 50 Hz sampling rate, which 
are post-processed by removing trend and outliers, most probably caused by instrumental 
error. Mean and standard deviation method was used to identify and remove data outliers, 
defined as all data points taking absolute value greater than three standard deviations about 
the mean value. Following the widely used procedure proposed by Van Dierendonck [13] for 
detrending signal amplitude and phase, the raw data were detrended with the high pass sixth- 
order Butterworth filter with the cutoff frequency at 0.1 Hz. As additional conditions, study 
contained only measurements made at elevation angles greater than 15°, to avoid errors 
induced by most intense multipath, and with continuous satellite contact or a time of lock 
greater than 180 seconds. Computed ionospheric indices used in tracking fluctuations of 
amplitude and phase are S,, defined as the standard deviation of the received signal amplitude 
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during 60 seconds of measurement, normalized to the average signal amplitude, and o,, 
defined as the standard deviation of the detrended carrier phase evaluated over 60 seconds. 


Hornsund Bay 


O : 
Presidente Prudente 


180°W 120°W 60°W 0° 60°E 120°E 180°W 


Figure 1. Geographical locations of receivers used in analysis — receiver set in high latitudes at Hornsund Bay, Spitsber- 
gen and in low latitudes at Presidente Prudente, Brazil. 


Two case studies on statistical analysis of the scintillating signals were undertaken during 
strongest geomagnetic storm in 2010, on 5-6" April, at Hornsund Bay, Spitsbergen, and 
moderate geomagnetic condition during 1* November 2011, at Presidente Prudente, Brazil. 
Geomagnetic disturbance index Kp is used as indicator of magnetospheric and ionospheric 
influence to the H component of geomagnetic field, which is in close relation with the gener- 
ation of ionospheric irregularities producing scintillations. Contributing parameters, as 
provisional 1 hour disturbance storm time (Dst) and 1 minute auroral electrojet (AE) indices 
are used in further data selection and in gaining detail information on all latitude magnetic 
activity. Data for all geomagnetic condition parameters were retrieved from the web archive 
of World Data Centre for Geomagnetism, Kyoto (wdc.kugi.kyoto-u.ac.jp). From Fig. 2 it is 
possible to notice that in case of 5-6" April 2010 level of Kp index reaching a maximum value 
of 8, while Dst and AE reached values of-81 and 2291 nT, respectively. While in case of low 
latitude data Kp index achieved not so high level with maximum value of 5, whereas Dst and 
AE come to-72 and 1261 nT, respectively. 


Analysis covered processing of 1, 3 and 12 hours data intervals for each of the tracked satellite 
PRNs, in order to test higher order moments on optimal amount of data mandatory for correct 
results. First analysis tests, not included in the paper, took data intervals of 3 and 12 hours, but 
due to averaging of large dataset, information from data points with pronounced scintillations 
had been lost and in all cases resulted to a close Gaussian PDF. Data analyses included 
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Figure 2. Fluctuations in geomagnetic Kp, Dst and AE indices during strong geomagnetic storm on 5-6*April 2010 
(right side graphs) and for moderate geomagnetic conditions on 1* November 2011 (left side graphs), showing initial 
and recovery phase of the storm. 


Gaussian, and Nakagami-m PDF of scintillating signals and its dependences on higher order 
moments. Nakagami-m distribution is known as mostly used in the statistical analysis for 
description of behaviour of strong ionospheric scintillation influence on the signal amplitude. 


Skewness is calculated as follows (third moment about the mean divided by sample variance 


powered by 1.5): 
aes) 
skewness = ——— (1) 
We 


and kurtosis as fourth moment about the mean divided by squared sample variance and 


normalized by 3: 


kurtosis= —— - (2) 
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n is number of values in a sample, x; represents i-th value in the sample and X is sample mean 
value. 


Choice of the PDFs in analyse is mainly relied on scientific and historical aspect [8, 14-20]. 
Nakagami-m distribution is known as mostly accepted solution in modelling of signal 
intensity, while Gaussian distribution is well known and widely used in statistical descriptive 
tool in science. In case of examination of signal intensity PDF, Gaussian distribution is used as 
reference distribution to measured data and Nakagami-m distribution. Gaussian distribution 
with zero mean, showed as acceptable, simple and effective solution for representation of the 
signal phase PDF for weak and moderate scintillations, is defined as 


P(I) = ROP ( e ) (3) 


where lis signal amplitude, o is a standard deviation and uis mean value. Gaussian distribution 
have found use across the various scientific fields, due to simplicity of use, good approximation 
of variety of natural phenomena and well clarified theory behind. 


Nakagami-m distribution, widely used in characterisation of ionospheric scintillation of signal 
intensity, radio links and wireless fading channels, is given by formula [9]: 


mz 
m my 2 


_ ml 4 
P= Tay expl- ae dl 


I represent single value of the signal amplitude calculated per 60 seconds, function shape 
parameter m is equal to 1/S,’, I(.) is Gamma function and <. > represents time average value 
of the signal amplitude. Spread parameter in Nakagami function, marked as <I>, is equal to 
signal amplitude mean value, which is set to 1. Nakagami distribution represents good 
approximation in describing multipath scattering with different groups of reflected wave [23]. 


Theoretical distribution of higher order moments, used in comparison with distribution of 
measured data, is calculated from the formula: 


KEPI "dI, (5) 


where P(I) is Nakagami PDF, I is signal amplitude and n=[1..4] is the order of the calculated 
statistical moment. From 3.381(4) formula of Gradshteyn and Ryzhik [23], integral was 
simplified to analytical problem. 


3. Results 


In order to obtain the information about probability distribution of scintillating signal, 
statistical analysis was performed by using higher order moments of signals phase and 
amplitude. Analysis was performed in two steps, first one included testing calculated higher 
order moments for received phase and amplitude for several time intervals (1, 3 and 12 hours), 
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and the second step included comparison of the PDF of measured data with PDFs obtained 
from theory. Initial test cases included data set lengths of 3 and 12 hours of measurements, 
fulfilling the condition that amount of data points should be large enough for correct calcula- 
tions of higher order moments. In these cases analysis showed good agreement with Gaussian 
distribution function, but only due to averaging over whole data set and losing valuable 
information about strong amplitude and phase fluctuations. This fact was confirmed with 
analysis of data samples of 1 hour length, which are shown in this paper. Additional condition 
is introduced, after detrending and removing outlier from measured data, only samples with 
more than 80% of usable data was taken in further processing. 
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Figure 3. Histograms of phase fluctuations made for one hour of observation, left and middle columns of graphs are 
for high latitudes measurements made on 5-6" April 2010, and right column of graphs for low latitude measurements 
made on 1* November 2011. Graphs in 2nd and 4th row represent PDF of measured data with fitted Gaussian distri- 
bution function. 
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Results of signal phase analysis for measurements in all three cases are shown in Fig. 3, with 
randomly chosen satellites PRN11, 4 and 2 (from left to right in figure) from group of cases 
with good alignment with Nakagami PDF in case of amplitude distribution, and PRN24, 29 
and 32 from group showing deviant characteristic. Most of analyzed samples showed good 
agreement with Gaussian distribution, but as it could be seen from the last row of plots in some 
cases (large skewness and/or kurtosis), phase distribution deviated from Gaussian PDF. 


Histogram of amprudo fluctuations - PRN11, 5th April 2010 Histogram of atmen fluctuations - PRN4, Sth April 2010 Histogram of amplitude fluctuations - PRN2, 1st November 2011 
2500 800 5000 
skewness = 0.50525 skewness = 0.62189 i 
E kurtosis = 9,3958 g kurtosis = 34311 2 praia 
3 2000 maen 54. 0,20568 7 3 maen S4 = 0.27025 3 4000 neces 
= 2 600 a maen $4 = 0.21683 
= E =] 
2 E] E] 
5 1500 5 5 3000- 
5 
E 3 400 El 
© $ D 
E 1000 £ E 2000 
S 5 3 
3 3 200 2 
E 500 E & 1000- 
5 3 3 
, i L - 
o 0 le mas. ba. 
0 0.5 1 15 2 o 25 3 4 
Signal amplitude [a.u.] nal a [au] Signal retik [au] 
2.5 15 r 2 
Data dist Data dist Data dist 
ol Gaussian |_ Gaussian Gaussian 
i Nakagami Nakagami 15 Nakagami |- 
3 A, 2 z = z a 
E] a 1 a 
315 i 5 3 
o ° o 
= j 2 24 
a 3 3 
S S 05 E 
a a a 
0.5 
05 
0 s o y 9 ! 
o 05 1 2 o 15 2 25 0 1 2 3 4 
Signal amplitude [a.u.] Signal amplitude [a.u.] Signal amplitude [a.u.] 
Histogram of one fluctuations - PRN24, 5th April 2010 Histogram of amplitude fluctuations - PRN29, 6th April 2010 Histogram dls fluctuations - PRN32, 1st November 2011 
4000 4000 T 2500 
= skewness = 0.97979 = skewnoss = 0.80501 i skewness = 1.4558 
€ kurtosis = 4.5568 = kurtosis = 4.741 = 2 kurtosis = 8.058 
2 3000 - mon S4 = 0.20677 a 3000 - maen S4 = 0.24612 3 ‘000 maen S4 = 0.25259 
E E E 
a S S 
5 5 5 1500- 
3 2000- 1 3 2000- 3 
a a a 
© o v 
E E E 1000 
3 3 5 
3 1000; 3 1000: 3 
£ £ & 500° 
5 5 5 
| IN i i ° I 
o -l AOS o al i T T 0 al Inn. Kanes 
0 25 0 05 1 15 2 25 0 3 4 
age edi [au] Signal amplitude [a.u.] Signal na [au] 
2.5 a 2.5 = 2 
Data dist Data dist Data dist 
ol Gaussian || ol Gaussian || Gaussian 
Nakagami Nakagami 15 Nakagami |- 
= z z 
2 2 ‘a 2 
3 15 3 15 3 
= = 21 
= = i = 
E E \ 1 g 
2 2 Ñ 2 
à o \ à 
\ 0.5 
05 05 N J 
0 ` ' ' =. o e ` ' 0 ' 
o 05 1 15 2 25 0 05 1 15 2 25 o 1 4 
Signal amplitude [a.u.] Signal amplitude [a.u.] Signal amplitude [a.u] 


Figure 4. Histograms of amplitude fluctuations made for one hour of observation, left and middle columns of graphs 
are for high latitudes measurements made on 5-6" April 2010, and right column of graphs for low latitude measure- 
ments made on 1% November 2011. Graphs in 2nd and 4th row represent PDF of measured data with fitted Gaussian 
and Nakagami distribution function. 


Results of analysis for amplitude scintillation are displayed in Fig. 4. Presented are PDFs of 
measured amplitude data with fitted theoretical Gaussian and Nakagami distributions. As in 
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Fig. 3 first two lines illustrates good agreement between experimental PDF and Nakagami 
fitted distribution. The agreement is more pronounced between 10:00 and 19:00 of local time, 
while the second case, showing deviations from Nakagami distribution, appears more in 
evening and early morning hours, especially in case of low latitudes. More detail analysis on 
larger data set is required for any further conclusions made on statistics of PDF deviation 
appearance and possible links to geo-physical parameters controlling these phenomena. 
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Figure 5. Representation of hourly mean higher order moments dependence of scintillation indices o,(first row) and 
S,(second row) for 5t, 6* April 2010 and 1% November 2011 watched from left to right side. Solid lines in the bottom 
graphs represent fitted theoretically derived skewness (blue) and kurtosis (red) from Nakagami distribution function 
for signal amplitude, which could be calculated only for values S, > 0.09. 


Fig. 5 depicts hourly values of higher order moments dependence on mean values of o, and 
S4. All the o, and S4 values are calculated every minute from randomly chosen PRN with full 
data set, which means that in specific hour PRN was visible all the time and measurements 
were performed without loose of lock. In graphs showing dependence from o,, skewness (blue 
y-axis and marks) and kurtosis (red y-axis and marks) are mainly taking values around zero, 
proving good agreement with Gaussian distribution. Only in case of 5" April 2010, it is possible 
to notice significant deviation from kurtosis, which occurs most possibly due to instrumental 
error. Bottom graphs displays higher order moments of amplitude, where skewness of 
measured data (circular markers) follows fitted distribution of moments (solid red line- 
skewness, and blue line-kurtosis) for smaller values of 54, while higher S4 skewness charac- 
teristic of measured data start rising much faster than theoretical one. Kurtosis controverts 
theoretical distribution, because of very sparse and chaotic behaviour, which could mean that 
Nakagami PDF is not ideal solution for representing amplitude distribution. 
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4. Conclusions 


Paper presents results of a GPS signal measurements statistical analysis at high and low 
latitudes. Focus was on the probability distribution of phase and amplitude under disturbed 
geomagnetic conditions for moderate and strong scintillations. Results show, as expected, a 
good agreement between the measured scintillating signal phase distribution and Gaussian 
distribution. More interesting case is for the signal amplitude, where a systematic increase of 
skewness and kurtosis with the S, values has been observed. In the case of weak scintillation 
(small S,) the probability distribution of amplitude is close to the Gaussian, while for strong 
scintillations (large S,) the skewness and kurtosis indicate considerable departure from the 
Nakagami distribution. This might be indicative of the non-Gaussian distribution of iono- 
spheric electron density fluctuations [19]. Dispersion in case of kurtosis could be due to the 
natural spread, but further experiments are required. 


Future research should consider more probability function models and x? goodness-of-fit tests 
for checking the precision of the distribution functions fit to the measured data. Also, more 
data sources, especially in low latitudes, and comparisons with in-situ measurements and 
other ground instruments measurements would be valuable reference in development of final 
model of scintillation signal distribution. 
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